
function [output,Loss_obsfit,Loss_seisfit] =rjMCMC(n_core,mcmcPar,seisPar,seisdata) 
%% 
storefile=strcat('result\Output_chain_test_',num2str(n_core));

firstrun=0;
%% ---------- parameters for rjmcmc updating --------------------
if firstrun==1
 
    % if any prior fractures revealed by downhole logs [x,y,strike]
    seisPar.fix_fracs = [
                          15	50	120.4;
                          15	63	28.6;
                        %  15	73.5	121;
                          115	37.5	120.2;
                        %  115	43.5	22.2;
                          115	58	120.3;
                          115	65	30.3;
                         % 115	71.5	114.4
                         ] ; 
                     
    seisPar.num_frac= randi([mcmcPar.kmin,mcmcPar.kmax]);   % initial arbitrary generated number of fractures
    seisPar.aperture = 10^(2.0*rand(1)-4.0);                % initial aperture for each fracture

    %% ---------------- initialization --------------------------
    n=seisPar.num_frac;
    [frac_info, s_error] = fracture_generation (seisdata, seisPar);
    fx=mcmcPar.obs/10;

    %% --------------- main iteration ---------------------------
    totaccept    = zeros ([mcmcPar.steps,1]); % acceptance ratio
    Loss_seisfit = zeros ([mcmcPar.steps,1]); % seismicity loss
    Loss_obsfit  = zeros ([mcmcPar.steps,1]); % tracer data loss

    birth_rec  = 0;
    death_rec  = 0;
    update_rec = 0;
    aBirth=0;
    aDeath=0;
    aUpdate=0;
    kt=0;
    tp=1;
else
    load(storefile)
    kt=size(output.frac,2);
    tp=round(kt*mcmcPar.thin);
    mcmcPar.steps = 200 ;                              % iteration times
%       totaccept    = vertcat(totaccept,zeros ([mcmcPar.steps,1])); % acceptance ratio
%       Loss_seisfit = vertcat(Loss_seisfit,zeros ([mcmcPar.steps,1])); % seismicity loss
%       Loss_obsfit  = vertcat(Loss_obsfit,zeros ([mcmcPar.steps,1])); % tracer data loss

end


for t=tp:tp+mcmcPar.steps-1
    t0=cputime;
% ----------------------------- move judgement -----------------------    
    Pr1 =exp(-s_error^2);
   % Pr1=abs(1/s_error);
    
    nbala= randi([mcmcPar.kmin,mcmcPar.kmax]);
    seisPar.num_frac=nbala;
    [~, s_error_new] = fracture_generation (seisdata, seisPar);
    if n < nbala
       Pr2=exp(-s_error_new^2);          % prior once adding a fracture  
       %Pr2=abs(1/s_error_new); 
    else
       Pr2=Pr1;                          % prior of more fracture
       Pr1=exp(-s_error_new^2);          % prior of less fracture
       %Pr1=abs(1/s_error_new);
    end
    birth=mcmcPar.arbc*min(1,Pr2/Pr1);
    death=mcmcPar.arbc*min(1,Pr1/Pr2);    
% ------------------------------------------------------------------   

% ---------------- fracture deletion/addition/updating -------------
    decision=rand(1);
    if decision <= birth && n < mcmcPar.kmax
        [ n,fracsPar,aBirth,seis_error, obs_error,fx] = frac_birth (mcmcPar, frac_info, s_error, seisdata, seisPar);
      %  fprintf('\t Birthrun ===============> %d \n',aBirth );
    elseif decision <= birth+death && n > mcmcPar.kmin
        [ n,fracsPar,aDeath,seis_error, obs_error,fx] = frac_death (mcmcPar, frac_info, s_error, seisdata, seisPar);
      % fprintf('\t Deathrun ===============> %d \n', aDeath);
    else 
        fracsPar=frac_info;
        n=size(fracsPar,1);
        aBirth=0;
        aDeath=0;
        seis_error= s_error;
      %  fprintf('\t Nothingrun ===============> %d \n ', 0 );
    end
    
    frac_info=fracsPar;
    s_error = seis_error;

    uLambda=rand(1);
    if uLambda>mcmcPar.lamda
      [ n,fracsPar,aUpdate, obs_error,fx] = frac_update(mcmcPar, frac_info);
      frac_info=fracsPar;
     % fprintf('\t Updaterun ===============> %d \n',aUpdate);
    end
    
    birth_rec  = birth_rec+aBirth;
    death_rec  = death_rec+aDeath;
    update_rec = update_rec+aUpdate;
% ------------------------------------------------------------------    

% -----------------------acceptance ratio --------------------------
    if t < 2
        totaccept(t)=(aBirth+aDeath+aUpdate)/t;
    else
        totaccept(t)=(totaccept(t-1)+aBirth+aDeath+aUpdate)/t;
    end
%------------------------------------------------------------------
    
% -------------------- store modeling fractures and fx ------------
    if mod(t,mcmcPar.thin)==0
        kt=kt+1;
        output.frac{kt}=frac_info;
        output.T{kt}=fx;
     %   fprintf('\t Storagerun ===============> %d \n', kt );
    end
    
    if mod(t,mcmcPar.update)==0
        save(storefile);
    end

%------------------------------------------------------------------
    Loss_seisfit(t)=seis_error;
    Loss_obsfit(t)=obs_error;
    
    fprintf('\t Iter = %d in %d core takes %10.3fs \n', t, n_core, cputime-t0);
end





 
  